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By elementary methods, I study the Love numbers of a homogeneous, incompress- 
ible, self-gravitating sphere characterized by a generalized Maxwell rheology, whose 
mechanical analogue is represented by a finite or infinite system of classical Maxwell 
elements disposed in parallel. Analytical, previously unknown forms of the complex 



o 

(N 

> 
O 

. shear modulus for the generalized Maxwell body are found by algebraic manipula 

> . 

tion, and studied in the particular case of systems of springs and dashpots whose 

strength follows a power-law distribution. We show that the sphere is asymptoti- 

^ ■ cally stable for any choice of the mechanical parameters that define the generalized 

Cf^ , Maxwell body and analytical forms of the Love numbers are always available for 

00 

I generalized bodies composed by less than five classical Maxwell bodies. For the 

homogeneous sphere, "real" Laplace inversion methods based on the Post-Widder 
. formula can be applied without performing a numerical discretization of the n~th 

derivative, which can be computed in a "closed-form" with the aid of the Faa di 

X 

' Bruno formula. 



I. INTRODUCTION 



Love numbers, named after A. E. H. Love 



23 



2J], represent a fundamental tool in geo- 



physics. From a physical standpoint, Love numbers basically represent properly normalized 
displacements and gravity potential variations in response to impulsive perturbations of a 
given harmonic degree. Since Love numbers for elastic Earth models can be easily gener- 
alized to the case of a linear viscoelastic rheology, they are useful to describe the response 
of the Earth on a broad spectrum of time-scales. As a consequence, using the Love num- 
bers technique, it is possible to address a number of relevant problems which range from 
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post-glacial deformations (see e. g., [4l| and re 
tions 45| , from post-seisrnic deformations 26 . 



'erences therein) to isostatic sea level varia- 

1 n 

31| to planetary tides and from Earth 



rotation instabi 
heterogeneities 



ities 19|, l28|] to the problem of dynamic compensation of internal mass 



35, 



42|. 



For an elastic, homogeneous, isotropic, incompressible and self-gravitating sphere, ex- 



tremely simple analytical forms exist for the Love numbers 1^, |28|], obtained from the 



solution of the Navier-Cauchy equilibrium equations 
displacement fields, and incremental gravity potential 



by. 



10 



an harmonic analysis of stress. 



22l |. The classic solutions pro- 



vided by Lame 201] and Thomson 48|| for the elastic compressible sphere and by Darwin 



for the viscous incompressible sphere have been later generalized to the viscoelastic, homo- 



eneous sphere 



29 



54 



55l | making use of the elastic-viscoelastic correspondence principle 



2l|. 



Amongst the existing closed-forms for the Love numbers, the one pertaining to the 
homogeneous Maxwell sphere has played a fundamental role during the past decades 29||, 
since the assumption a Maxwell viscoelastic rheology largely explains some of the geo- 
)hysical observations accompanying post-glacial rebound and long-term mantle dynamics 



30l |. Current investigations in the field of global geodynamics, however, are performed 



using multi-stratified Earth models compatible with seismological evidence, in w 



lich the 



equilibrium equations are generally solved assuming a complex viscosity profile [46j , whose 
depth-dependence is varied until surface observations (geodetically observed deformations, 
relative sea level and gravity field variations) are satisfactorily reproduced 30 1. 

In this work, we go back to the homogeneous, incompressible and self-gravitating, vis- 
coelastic sphere (hereinafter if-sphere), to discuss some aspects that have been apparently 
unnoticed so far, possibly because of the large success of the simple (but simultaneously 
realistic) Maxwell rheology, and of the ensuing numerical applications to multi-layered 
models. In particular, we extend the Love numbers formalism to the case of generalized 
(discrete) Maxwell bodies (hereinafter GMBs), whose properties are of particular interest 



in various fields of physics 



. In general, a GMB results from 



25l | and geophysics 

the one-dimensional arrangement of various classical Maxwell bodies (CMBs), whose ma- 
terial parameters are chosen so that to reproduce physical (or geophysical) observations 
[sl. Here we limit our attention to discrete finite or infinite GMBs obtained by elementary 
parallel arrangement of CMBs, and we address the problem of Laplace inversion of the 
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so-generalized Love numbers. More complex combinations of CMBs, such as the ladder 



networks, provide fractional constitutive re^ 



ationships which are of particular interest 



in the theory of electromagnetic systems [iJ] . Love numbers spectra corresponding to these 
arrangements will be considered elsewhere, in view of possible geophysical applications. 

The paper is organized into four sections. After reviewing in Section Hi Al the properties 
of discrete GMBs composed by a finite number of CMBs, we consider the complex shear 
modulus of a discrete, infinite GMB with mechanical parameters distributed according to 
a power-law, also giving - apparently for the first time - closed forms for the viscoelastic 
material functions in terms of classic special functions, reported in Section III B[ Then, in 
Section fll CI the Love numbers for the if-sphere are generalized to a rheology described by 
finite GMBs, also discussing their Laplace-inversion by means of traditional methods. In 
the final part (Section III Dl) . we address the problem of Laplace inversion of the generalized 



Love numbers by means of Post's formula 3J]. Seen the simple structure of Love numbers 



in the Laplace domain (this is a consequence of the geometrical simplicity of the if-sphere), 
"closed forms" are available for the second-order Bell polynomials that enter the Faa di 



Bruno formula 



171], hence, in principle, the n-derivative of the Love numbers - required 



in Post's formula - is available analytically. 

II. RESULTS 
A. Discrete GMBs 

The classical Maxwell body (CMB) is a simple mechanical system composed by a spring 



connected in series with a dashpot 



25 



33| . The quasi-static creep or relaxation of the 



CMB can be studied in the Laplace-transformed domain introducing the complex shear 
modulus 

where s = x + iy is the complex Laplace variable and the material parameters /i (/i > 0) 
and rj [rj > 0) represent the rigidity and the viscosity of the spring and of the dashpot, 
respectively. The ratio 

r = ^ (2) 

is Maxwell relaxation time of the CMB. 
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Function /i(s) fully describes the response of GMB, expressed by the stress-strain rela- 
tionship [5]. The creep compliance J{s) and relaxation modulus G{s), which represent the 
response of the GMB to a unit stress and strain, respectively, are in fact related to fi{s) 
by G{s) = 2jl{s)/s and J{s) = \/2sjx{s) (see e. g. 25|). Functions J{s) and G{s), also 
referred to as material functions of the GMB, are not independent one from each other. 



being linked by the reciprocity relation J(s)G(s) = (^S- |25l|). 

By the combination rule for mechanical analogues j^, [s^ , the complex shear modulus 
of a discrete GMB composed by N GMBs disposed in parallel is 

N 



Ks) = ^lln{s) (3) 



n=l 

where, from ([1]), the complex shear modulus of the of the n— th GMB is 
with rigidity fin > ^ and viscosity rin > 0. The constant 

Tn = (5) 

represents the Maxwell relaxation time of the n-th GMB component (hereinafter, it will 
be assumed that times r„'s are distinct). 

In terms of r^, the complex shear modulus of a A^-elements GMB reads 

N 

n=l ^/^n 

showing that /i(0) = and that fl{s) has exactly isolated poles for s G M", located at 
Sn = -l/^-n- From 

djl{s) ^ fln/Tn 



E H'n/ In /-N 



and by the positivity of fin and it follows that fi{s) is strictly monotonic for s G 
These properties show that the N zeros of fl{s) are interlacing the poles in s G Mg 13 
Since the A;-th derivative of the complex shear modulus is 

N I 



/x(s) is a (7^ function for s G (i. e., it is infinitely different iable along the real positive 

QviG^ ^irViir»Vi tinanytia 1"Viti QTM^lir>QV\ili1""\7' r\f +Viti "T*tiQl" Pr\a+ — \hJ\AApnv T,QT^lQr»ti imrtiyair^n Tnti1"Vir\r] 
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to the Love numbers problem for the homogeneous sphere, as we will discuss in Section 
III Dl below. In addition, since 

{-If fi'^'\s)<0, seR+, (9) 



we note that fi{s) is a completely monotonic function (e. g. |25|). 

The limit of for c)0 is not straightforward. For instance, it is clear that an 

infinite GMB composed of identical springs = /iq) and dashpots (?7„ = r]o) combined 
in parallel does not have a finite complex shear modulus (i. e., series is divergent). 
This shows that finite values of /i(s) can be obtained only with appropriate combinations 
of elastic and viscous elements, with varying strengths. A case study will be investigated 
in the next section. 



B. A power law, discrete GMB 

We consider, study, the response of a GMB with moduli following a power-law 

distribution, with 

f^n = ^, peN, fi*> 0, (10) 
Tin = -. g e N, r/* > 0, (11) 

where /z* and rj* are a reference rigidity and viscosity, whose ratio defines the time constant 

r* = ^. (12) 

The two-parameters GMBs described by (11 01) and (II II) are particularly useful since 
closed-form expression are available for the complex shear modulus in the case = oo, as 
we will show below. This implies, in particular, a closed-form for the material functions 
J{s) and G{s)^ which are generally not available for finite arrangements of mechanical 
analogues. At the same time, a power-law distribution of material parameters is sufficiently 
general to be potentially useful for numerical applications in physics and geophysics. An 



example has been recently given by Spada 



461 ]. who has employed this distribution to 



study the Love numbers of a multi-stratified Earth model and has anticipated one of these 
analytical forms in the particular case (p = 0, g = 2). Here the mathematical aspects 
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presented in 46(] are considered more in detail and extended to any value of the integer 
exponents p and q. The case {p, g) G M will be investigated in a follow-up study. 
It is now convenient to normalize the complex shear modulus 

m{s) ^ (13) 

which, using (fTU]) and (fTTj) with and (jl]), gives 

N 

m(z;p,q) = 'Y^ , (14) 

with 



nPz + ni ' 

n=l 



z = sr* e C. (15) 



For finite values of and arbitrary distribution of moduli, the series (iMl) cannot be 
summed to provide a closed-form complex shear modulus. However, a general result that 
can be easily established valid for all N values (including N = oo), is 

m{-;q,p) = -m{z;p,q), (16) 
z z 

showing that the modulus of a given GMB can be obtained from that of a complementary 
GMB, in which springs (with distribution determined by p) and dashpots (g) are inter- 
changed. As a consequence of the symmetry-duality relationship flTB]) . the summation of 
f lT^ can be limited to p < q. 

For a GMB composed by an infinite number of CMBs, the normalized complex shear 
modulus is 

M{z;p,q) = lim m{z;p,q), (17) 
with m{z;p, q) given by ( JT^ . Hence we are interested in the study of the series 

oo 

M(z;p,g) = V ^ P,qeN, (18) 

ra=l 

for which the conditions of convergence (divergence) are the same as for the series 
^^-^ 1/ (n^z + ra'?) . Since l/ln^z + n'^\ < and Yl'^=i ^ / ''^'^ convergent for q > 2, 

by the Weierstrass M-test for the series of complex functions (see e. g. ISj), the (uniform) 
convergence of ( fT8l) in this range of q values is proved. By a similar argument, it can be 
easily shown that a further condition of convergence is p > 2. Hence, we conclude that 
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sufficient condition for the uniform convergence in the whole complex plane of M[z;p,q) 



IS 



{p,q)eA, A = {p>2}U{g>2}. (19) 

The divergence of (1181) for {p, q) ^ A can be shown in a straightforward way. 
The poles of M{z;p, q) are found at 

Zn = -n''-". (20) 

hence they are simple and, for p ^ q they are countably infinite (in the particular case 
p = q, the infinite GMB degenerates into a CMB with Maxwell time r*, with M{z]p,q) 
showing a single pole Zi = —1). For any p and q value, the poles are Zn G M^, and, from the 
general properties of complex modulus /i(s), discussed in Section Hi At they are interlaced 
with the zeros of M{z; p, q). Points z = —oo and z = are accumulation points of poles for 
q > p and q < p, respectively. It is also of interest to observe that, in the limit for z oo, 
M{z] p, q) is only determined by the strength of the springs (this is physically sound, since 
the limit z ^ oo corresponds to the small times limit). In fact, from (ITSl) one obtains 

lim M(z;p,g) = C(p), (21) 



where C, is Riemann zeta function [l|. Hence, M{z;p,q) is bounded at ^ = oo only for 

p e A. 

With the help of tables of series [l^ and of an algebraic manipulator, it is straightforward 
to verify that closed-form expressions for M{z; p, q) exist in the case of discrete GMB with 
N = oo. As discussed in Section III Al they can be used to obtain closed-forms for the 
material functions J{s) and G{s) of the GMB, which are usually not available for finite 
values of N. These analytical formulas are useful since they allow for a compact expression 
of M{z; p, q) but their complexity, also manifest from the infinite number of poles and 
the presence of accumulation points of poles along the real negative axis, can make the 
Laplace inversion of Love numbers in the time domain practically problematic, as it will 
be discussed in Section Hi CI 

The closed-form expressions that can be obtained by Equation f|T8|) involve classical 
special functions (the derivative of the digamma function il){k, z) and the Riemann zeta 
function C,{s), respectively), as illustrated in Tabled for low values of p and q. Definitions 
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and elementary properties of these functions are found in e. g. Compact forms of 
M{z;p, q) however exist also for larger values of p and q. For instance, have verified that 

with il>{k^ z) = -^ijj{z) where ip{z) is the digamma function and $,kik = 1, . . . q) are solutions 
of the algebraic equation z + [1 + = 0. By virtue of the reciprocity relationship f|T6l) . 
the complex modulus M{z; q, 0) {q > 3) can be easily determined from (l22!l . 

Though we have only studied function M{z] p, q) for a limited number of p and q values, 
we conjecture that algebraic manipulation can provide 'closed-forms" for any value of 
parameters p and g, though these formulas could be too complex (and the CPU time 
required for manipulation exceedingly long) for being of any practical use. 

TABLE I: Normalized complex modulus M(z;p,q) (jlSp for some GMBs wit power-law distri- 
bution of moduli, characterized by low values of p and q, with {p,q) G A (jl9p . Here 7 is Euler 
constant (7 = 0.577215...), '4){k,z) is the fc-th derivative of the digamma function, defined as 
= V {z) /T{z) where T{z) is the Euler gamma function, and C,{s) = X^fcLi ^ is Riemann zeta 
function p|. Note that M{z;p, q) obeys the reciprocity-symmetry relationship (I16p . valid for any 
N. z = Zn gives the location of the poles of M{z;p,q). For (p, q) = (0,2), M{z;p,q) can be 
obtained from published tables (see e. g. Equation 4. /1. 421 of |l3l]), while for p = q, it follows 
from the definition of Riemann function [l|. All the other forms have been obtained by algebraic 
manipulation. 



p 


q 


M{z;p,q) 


Zn 





2 


V2(— 1 + TT^/z cothvr^/i) 
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72 (-2 + VTA/icoth — ) 
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z(7 + V(0,l + -)) 

z 


1 

n 


p 


p 


zC{p) 
l + z 


1 



9 



C. Generalized Love numbers for the if sphere 

At a given harmonic degree £, the Laplace-transformed Love numbers for the if-sphere 
can be cast in the form 



where Lf{t) is the "fluid" limit of the Love number (i. e., Lf{t) = limg^Q L{s)) , p,{s) is 
the complex shear modulus of the CMB (or more generally, of the GMB) that mimics 
the rheological behavior of the sphere and fie is the elastic rigidity of the sphere. With 
appropriate functions Lf = Lf{i), Equation (1231) is useful to describe vertical and horizontal 
component of displacement, and the incremental gravitational potential, for Love numbers 



19|, 



28|, 



4l|. In Equation fl23l) . I have introduced the non- 



of both tidal and loading type 
dimensional constant 

A= = ^Jl±£±l^, (24) 

i pga 

where p is the density of the sphere, a is its radius, and g is gravity at the surface {g = 
y^nGpa, G being Newton gravity constant). At a given degree i, is a measure of the 
ratio between elastic stress (governed by u^) and gravitational stress (described by pga). 



For the "average" Earth, p^l PQO. ~ 0.60 [57 1. 



By substitution of (jlj) into (l23l) . the Love numbers for a //-sphere with a GMB rheology 
can be easily studied. The poles of L{s), which correspond to the zeros of 1 + X^fi{s)/pe, 
are all real and negative. In fact, recalling from Section III Al that p{s) is monotonic for 
s G M, vanishes for s = and has — 1 more zeros for z G M~, the zeros of 1 + \^jl{s)/ pe 
must be found for s G M^, being shifted to the left relative to those of p{s) because of the 
additive term " 1" . 

From above, we conclude that any (incompressible) if-sphere with GMB rheology is 
stable with respect to surface or tidal loading, for perturbations of any harmonic degree 
and regardless of its material properties. This also holds for N = oo, since adding more 
CMBs to the system would not change qualitatively the distribution of the zeros of p{s). 
This stability property is certainly violated for compressible spheres of initially constant 



density, as clearly illustrated by 15| in the case of a simple CMB. 

For a GMB composed of A^ elements, an analytical Laplace inversion of L{s) can only 
be obtained, in nrincinle. for A" < 4. This can be seen bv substitution of fllll into (12311 . 
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which provides 



F(s) = 



Lis) 
Lf 



N 



(25) 



^ s + 1 



where fi'^ = fin/f^e- Hence 



where 



F(s] 



Pis) 

Q{s) 



Pis) = fl(s + ^) 



and 



AT 



Q{s)=Pis)+\'J2^''nS 



n=l 



^ ( 1 

n 



are degree N polynomials in the variable s. 

Hence, by the Heaviside expansion theorem (see e. g., 
number for the GMB can be cast in the multi-exponential form 

N 

Lit)=LJit) + Y,Lne'-\ t>0, 



(26) 



(27) 



(2^ 



the time-domain Love 



(29) 



n=l 



where 6 is Dirac's delta, s„ [n = 1, . . . N) are the (real and negative) distinct roots of the 
algebraic equation 

Q{s) = 0, (30) 



and elastic and viscoelastic components of Love number are 

Pis) 



lim 

s^oo Qis) 



and 



PiSn 
Q'iSr. 



n = l,...N, 



(31) 



(32) 



respectively, where Q' is the first derivative of Q{s). An exact solution of Equation (!30|) 
is only possible analytically for < 4, since by the Abel-Ruffini "irnpossibility theorem", 
general quintic equation cannot be solved in terms of radicals (e. g., 47|). We remark that 
times —1/sn bear no obvious relationship with the time constants r„, defined by ([5]). 

The existence of closed forms for the Love numbers for GMBs with A^ < 4 guaran- 
tees the possibility of obtaining analytical results for particular GMBs of great interest in 
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geophysics. This is the case of Burgers rheology, a four-parameters model which is tradi- 
;ionally represented by a CMB combined in series with a Kelvin- Voi gt e lement (see e. g., 
33|), and widely employed in the study of post-seismic deformations [26|, l31|, post-glacial 



rebound 



18 



27 



37|, 



561 ] and planetary dynamics [53]. Since it has been shown that such 



disposition is mechanically equivalent to a four-elements GMB composed of two CMBs in 
parallel , a closed-form expression of the type (!29|) with = 2 is certainly possible for 
the Burgers if-sphere, where the explicit relationship between L^, Ln and (z = 1, . . . A^) 
and the four free parameters of the Burgers body (/ii, /i2, ?7i, V2) can be obtained by lengthy 
algebra, since a quartic equation is involved. 

For > 5, the Laplace inversion of the Love numbers can be only performed by a 
numerical evaluation of the roots of polynomial Q{s) in Equation (l28l) . again followed by 
the application of Heaviside expansion theorem. The multi-exponential form given by fl29|) 
is therefore still formally valid for > 5, but the coefficients cannot be expressed explicitly 
in terms of the mechanical parameters of the GMB. 



D. Post Widder formula and the sphere 



The simple analytical structure of the generalized Love numbers ( l25l) allows, at least 



formally, alternative approaches to the 



the Post-Widder (PW) formula [32 



51 



ap 



ace inversion, based on "real" methods such as 



521] (a modern, detailed proof of Post's inversion 



formula can be found in l3|, with a nice comment on the ill-posedness). In numerical 

n 

applications (e. g., [49;]), the main advantage of PW formula is that it does not require 



root-finding numerical algorithms, which can become unreliable especially for large A^, 
when equation ( l30l) may possess densely packed (and thus numerically difficult to resolve) 
roots on the real negative axis [4^. In the context of this study, as we have discussed in 
Section riIB| for A^ 1-^ 00, the roots are countably infinite and accumulation points of poles 
appear, that enhances the numerical difficulties. 

The PW formula requires the computation the derivatives L^"'\s) along the real positive 
axis (hence the attribute real) and the evaluation of the limit of a sequence according to 

,(,.U.ti)!f?y'"i<..f;) (33, 



TO— >oo n\ 



51 



52| requires L(s) G C°° for s G M+ 6|. The convergence of sequence (j33il is logarith- 
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mically slow, but it can be efficiently accelerated [49|, |50| without seriously compromising 
the performance of numerical computations - at least in the geophysical applications per- 
formed so far 44, [4^. Lacking, in general, an analytical expression for L^^\s), numerical 
application of the PW formula requires a finite-difference discretization, a noisy numerical 
operation that demands a multi-precision environment (a nice tool is offered by FMLIB 



40( 1) to prevent the phenomenon of catastrophic cancellation [39]. As we will show below. 



can be evaluated analytically in the present context, thus avoiding numerical the 
discretization which constitutes a major limitation of the PW method. 

Application of the PW inversion method to the Love number problem for the if-sphere 
is feasible, since L{s) is smooth (i. e., L{s) G C°°) for s G M'*', being fi{s) itself smooth in 
this interval. Writing 

F(s) ^ ^ (34) 



gives 



with 



F{s) 



g = g(s) = A 



Lfii) 
1 

2/i(s) 



(35) 



(36) 



The n-th derivative of F (5), required in Equation (!33|) . can be expressed using the Faa di 



Bruno chain rule formula 
Namely 



17 



36l | for the derivative of the composite function F = F{g{s)). 



FW(.) = 5^FW(^)i?„,,(^«,^(2),...,^ 



(n-k+l)\ 



(37) 



k=0 



where, using (135!) . the k-th derivative of F with respect to g is 



F^'\g) = (-l)^T^ 



+ 



\k+l 



(38) 



and Bn,k denotes the incomplete Bell polynomials (also known as second kind Bell poly- 
nomials) [2I, |7| , defined as 



where the sum is over all sequences of non-negative integers ji, j2, ■ ■ ■ jn-k+i which are 
solutions of equations ji + j2 + ■ ■ ■ = k and ji + 2j2 + Sj^ + . . . = n. 
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In the present context, a "closed-form" expression for Bell polynomials can be obtained 
recalling that for a GMB with elements, the m-th derivative of fl{s) is given by Equation 
(IH]). Hence, for any integer m, 

y V ; \ ) ^As + l/r„ "^+1 ^ ' 

11=1 ^ ' ' 

can be used in the right hand side of Equation (1371) . obtaining a fully explicit (but extremely 
complex) expression for F("^(s). In this way, the major shortcoming of the PW formula, 
namely the numerical noise amplification produced by repeated differentiation of I/(s), can 
be circumvented for the if-sphere (but to the cost of very complex algebraic computations). 
Therefore, the generalized Love numbers for the iJ-sphere can be expressed as a limit of the 
sequence ( 133|) . which in principle may constitute an alternative to the classical root-finding 
approach, especially for large values. 



III. CONCLUSIONS 



Our main conclusions can be summarized as follows, In the case of discrete GMBs 
composed of an infinite number of CMBs disposed in parallel, analytical forms for the com- 
plex shear modulus are available in the case of material parameters distributed according 
to an (integer) power-law (see Equation [TOj) . These forms involve classic special functions, 
and are moderately simple for low values of the powers. After algebraic manipulation of 
several case studies, we conjecture that analytical (but exceedingly complex) moduli can al- 
ways be formally determined, ii) For finite GMBs composed by limited number of elements 
(in particular, A^ < 4), the Love numbers of the if-sphere can be determined in closed 
form. These Love numbers are asymptotically stable for any value of A^, provided that the 
iJ-sphere is incompressible. For A^ > 5, standard numerical instruments can be used to 
determine the poles of the Love numbers, which could however suffer from the presence of 
accumulation points for the poles. Numerical difficulties in the numerical Laplace inversion 
of the Love numbers are well documented even in the case such singularities do not enter 
into play 4^. in) The extremely simple algebraic form of Love numbers for the if-sphere 
allows for a closed-form construction of the Bell poljTiomials, which enter the Faa di Bruno 



formula for the n-th derivative [17|. Therefore, the numerical difficulties that follow from 
the numerinal discretization nf the derivative ll J. Isol . can be nartlv circumvented. 
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